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“Alchemical” interpolation paths, i.e. coupling systems along fictitious paths that without realistic 
correspondence, are frequently used within materials and molecular modeling and simulation pro¬ 
tocols for the estimation of relative changes in state functions such as free energies. We discuss 
alchemical changes in the context of quantum chemistry, and present illustrative numerical results 
for the changes of HOMO eigenvalues of the He atom due to a linear alchemical teleportation- the 
simultaneous annihilation and creation of nuclear charges at different locations. To demonstrate 
the predictive power of alchemical first order derivatives (Hellmann-Feynman) the covalent bond 
potential of hydrogen fluoride and hydrogen chloride is investigated, as well as the van-der-Waals 
binding in the water-water and water-lrydrogen fluoride dimer, respectively. Based on converged 
electron densities for one configuration, the versatility of alchemical derivatives is exemplified for 
the screening of entire binding potentials with reasonable accuracy. Finally, we discuss constraints 
for the identification of non-linear coupling potentials for which the energy’s Hellmann-Feynman 
derivative will yield accurate predictions. 


INTRODUCTION 


Ever since the introduction of Hess’ law and Carnot’s 
cycle, chemists have known that some properties, called 
state functions, always change by the same amount when 
a system is moved reversibly from one state to another— 
regardless of how the change has been implemented. The 
freedom to choose any paths, even paths without any 
realistic correspondence except for the endpoints, is ex¬ 
ploited within many applications. We generally refer to 
“alchemical” paths as paths that cannot be followed and 
verified through experimental observations. For exam¬ 
ple, Fig. [TJa) illustrates how, according to Hess’ law, the 
change of enthalpy of reaction can be calculated either 
by following the (realistic) reaction path, or, just as well, 
by following a more convenient yet non-realistic (alchem¬ 
ical) reaction path to product via dissembled elemental 
states as intermediates. Depending on the choice of state 
function, external conditions, system, and process, real¬ 
istic reaction paths can be significantly more challenging 
because they can involve many intermediate and transi¬ 
tion states which are difficult to identify and characterize. 
Even worse, they might even be experimentally impos¬ 
sible to probe, as it is the case for the chemistry of the 
earth’s core, some other planet’s bio-sphere, for distant 
historical or future events, or for very slow or very fast 
processes. 


Within the atomistic theories of quantum and statis¬ 
tical mechanics, any path connecting the Hamiltonian of 
some initial molecule or material system, Hi, to some 
final system Hf, can be defined in a coupling order pa- 



FIG. 1. Alchemical cartoons, (a) The same enthalpy change, 
Ah, is obtained for a realistic (full) or an alchemical (dashed) 
coupling between initial (i) and final (/) states as a func¬ 
tion of reaction progress £. (b) Alchemical paths connecting 
compounds on two two-dimensional binding potential energy 
surfaces corresponding to two different stoichiometries, {Zi} 
(white) and {Zjj (gray), respectively. Having calculated Et, 
for some initial system i (filled circle), alchemical paths couple 
to energies (open circle) of different geometries / with same 
stoichiometry (dotted), different stoichiometries with same 
geometry /' (full), or different geometries and stoichiometries 
f" (dashed). 

rameter A as long as as the end-points are met, HHHji.e. 
(Hi, A = 0, 

H( A) = J H x , 0 < A < 1, (1) 

i H f , A = 1, 

where 0 < A < 1. H x in Eq. 0 denotes some intermedi- 
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ate state at A, not necessarily differentiable. At bound¬ 
aries of first order phase transitions, for example, the 
entropy (state function) is not continuous in tempera¬ 
ture (A). Often, H( A) is (arbitrarily) chosen to be lin¬ 
ear in A, i.e. H( A) = Hi + \{Hf — Hi). As alluded to 
above, H\ does not have to be realistic for all values of 
A. Thermodynamics textbook examples of such changes 
include the calculation of the errors made when relying 
on the ideal gas equation. Introduced as “computational 
alchemy” mm in the realm of computational chemistry, 
this concept has successfully been used for the interpola¬ 
tion of forces and energies for molecular dynamics (MD) 
and Monte Carlo (MC) simulations. Also for the pur¬ 
pose of quantum mechanical observables, we can denote 
any such unrealistic path as “alchemical”. mm We note 
however that Eq. (|T|) is also known as “mutation path” 
or “adiabatic connection”. CHS] 

An even more intriguing possibility for exploiting the 
freedom of alchemical changes relates to the challenge 
of rational compound design (RCD). RCD attempts 
to circumvent (or at least reduce) the combinatorially 
scaling challenge of having to virtually enumerate and 
screen larger subsections of chemical or materials com¬ 
pound space using computationally demanding simula¬ 
tion methods. It has already been shown to yield promis¬ 
ing results for the virtual atomistic control of material, 
nanoparticle, and even molecular structures. (It), [TT] Be¬ 
cause of the vastness of chemical compound space (CCS), 
identification of novel compounds that meet desired prop¬ 
erty requirements still remains a challenge. Once an 
alchemical interpolating path, H{ A), is defined, property 
derivatives with respect to A can be evaluated m (see 
Sec. [f. Similar to an iterative gradient descent-like al¬ 
gorithm, one can thus navigate gigantic combinatorial 
compound libraries at dramatically reduced computa¬ 
tional costs by visiting the most promising compounds 
one after the other while avoiding the least promising 
candidates. [HI ITS] 

The concept of connecting different systems via Eq. ([l]) 
has been in frequent use in various research fields, in¬ 
cluding computational engineering, physics, biophysics, 
and chemistry. Here, we first briefly summarize the most 
common application of Eq. ([I]) to calculate free energy 
changes, or alloy formation energies in Sec. [] In Sec. [] 
we review the quantum mechanical treatment of alchem¬ 
ical changes. To this end, we mainly rely on the use of 
density functional theory (DFT) even though analogous 
arguments can be made using conventional wave-function 
based quantum chemistry methods. In Sec. j] we present 
numerical results that demonstrate the use of alchemical 
derivatives for the screening of entire potential energy 
binding surfaces with semi-quantitative accuracy with¬ 
out additional self-consistent held calculations. 


COMMON ALCHEMICAL APPLICATIONS 

Free energy is one of the most important state func¬ 
tions in chemistry. Since it is a statistical average, large 
numbers of configurations need to be taken into account 
to yield accurate predictions. m E.g., calculating a free 
energy of solvation following a path that mimics the real¬ 
istic complex process of reversible microscopic immersion 
of the solute into a condensed ensemble of a very large 
number of solvent molecules would imply a severe simu¬ 
lation effort that ensures that all relevant degrees of free¬ 
dom have sufficiently been sampled. Furthermore, to ac¬ 
count for hysteresis effects, this simulation should be re¬ 
peated for various initial conditions and immersion rates. 
And one would have to start anew for any changes made 
to temperature, pressure, or solvent and solute species. 
Alternatively, one could also calculate the change in free 
energy with respect to some solute for which the free 
energy of solvation is already known. Thermodynamic 
integration, i.e. numerical integration of the statistical 
mechanical average of the “alchemical force” along the 
path converting known solute (A = 0) into query solute 
(A = l),i 

^-jCX^irX- (2) 

Jorgensen and Ravimohan|7] proposed an even more 
efficient alternative: One can also estimate the change 
in free energy of solvation due to changing the solute 
using perturbation theory and MC simulation. Specif¬ 
ically, they considered the effect on the free energy of 
hydration due to an alchemical change of a methyl into 
hydroxy-group, AG = G f - Gi = G C h 3 ch 3 - G C h 3 oh- 
One can show that if the sampling of the two states, Hi 
and Hf , yields sufficient overlap, the corresponding free 
energy difference can be accurately predicted using per¬ 
turbation theory, 

e-P AG « . (3) 

Here, 1//3 = ksT, and the right-hand-side refers to 
the average of the Hamiltonian difference Boltzmann’s 
weight over a trajectory generated using Hi. The authors 
used a linear interpolation of force held parameters for 
methanol and ethane, H( A) = JJch 3 oh — A(Hch 3 ch 3 — 
IJch 3 oh), from which the energy can be calculated for 
any A. 

As such, alchemical changes enable the prediction of 
changes in free energy differences without having to ac¬ 
tually model the realistic process under investigation. 
Linear interpolation approaches have been applied to 
free energy calculations in various chemical and biologi¬ 
cal systems. TOD] Smith and van Gunsteren found that 
non-linear alchemical coupling not necessarily leads to 
linear free energy changes. ED Further applications of al- 
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chemical coupling to the estimation of free energy differ¬ 
ence include the free energy of hydration of ions using ab 
initio molecular dynamics, [22] differences in free energy 
of binding between various host-guest complexes,[53] free 
energy differences at phase boundaries to predict melting 
points, P3J [55] the free energy of mixing to identify eutec¬ 
tics in ternary mixtures of molten alkali-nitrate salts. [5H] 
kinetic isotope effects,[27] as well as constraints on the 
composition of the Earth’s core. [53 

But also from the solid state point of view the con¬ 
cept of alchemical coupling is used for the prediction of 
properties of disordered materials, such as co-crystals, 
solid solutions, or solid mixtures, as a function of mole- 
fraction. [29] It is computationally difficult to deal with 
such mixed disordered systems since the minimal self¬ 
repeating units can become very large. As a result it 
is nearly impossible to set up disordered systems within 
periodic boundary conditions. One alternative consists 
of using cluster-expansion methods [30], another alter¬ 
native, akin to alchemical coupling, is the virtual crys¬ 
tal approximation (VCA) [31] which averages the system, 
rather than explicitly representing the full system. One 
of the simplest disordered class of materials are ternary 
semiconductors, A^Bi-^C, where AC and BC are two 
different semiconductors while x is the mole-fraction be¬ 
tween A and C. Consider, for example.[35] Eq. 0 applied 
to Al^Gai—g,As. H(x ) — a s -\-x(Haw s RGaAs)- The 
linear interpolated alchemical path describes an aver¬ 
aged Hamiltonian between AlAs and GaAs for any mole- 
fraction of A1 and Ga. 

ALCHEMY IN QUANTUM MECHANICS 
Fictitious systems 

Within a first principles notion of CCS, m one can 
view every compound in any geometry as a state de¬ 
scribed by a unique Hamiltonian H. More specifi¬ 
cally, the total potential energy’s molecular Hamilto¬ 
nian, H, is a function of a given set of nuclear coordi¬ 
nates, charges, and number of electrons, {R i,Zi,N e }, 
respectively. Without any loss of generality, we here 
rely on the Born-Oppenheimer approximation, neglect¬ 
ing all non-adiabatic electronic or nuclear quantum ef¬ 
fects. Studies of alchemical paths have historically pro¬ 
vided essential insight into the density functional the¬ 
ory (DFT) formulation of the many-electron problem in 
molecules. [331 kSJ In 1974, Harris and Jones introduced 
an adiabatic connection, [35] coupling the system of inter¬ 
est to an fictitious but relevant system of non-interacting 
electrons, 

H(X) =T + XV ee + V ext , (4) 

where T, V ee , and V ex t represent kinetic energy, electron- 
electron interaction energy, and external potential en¬ 


ergy operator. By changing A from 1 to 0, one can 
dial in the electron-electron interaction. For A = 0, 
the electronic Schrodinger equation can thus be solved 
analytically, providing useful information on properties 
such as the exchange-correlation hole. ]36fi351 an impor¬ 
tant ingredient for current exchange-correlation poten¬ 
tial development efforts. (3HIH5] Another important study 
of electron-electron interaction, carried out by Seidl, 
Perdew and Levy, introduces the limit of strictly cor¬ 
related electrons. @3] Replacing the variable A = ^ for 
0 < n < 1 in Eq. Q, one obtains a coupled system 
where electron-electron interaction is dominant. 

E. B. Wilson introduced the idea to alcliemically cou¬ 
ple any system to the uniform electron gas. Based on 
this path, he derived an expression for an exact four¬ 
dimensional density functional theory, integrating over 
three spatial and one A-dimension. [33 m] Subsequently, 
Politzer and Parr[45] showed that, by defining free-atom 
screening functions, Wilson’s functional can be decom¬ 
posed into kinetic and potential energy of N e electrons. 
These definitions of DFT related alchemical paths con¬ 
stitute the underlying framework for the results and dis¬ 
cussions here within. 

Within DFT, [33] we can explicitly calculate E( A) for 
any iso-electronic change of geometry and composition, 
i.e. under the constraint that f dr n x (r) = N e V 0 < A < 
1 , 

E[n x ,X] =T[n x ] + V ee [n x ] + j dr n x (r) v ext (r, A). (5) 

Here, the coupling is introduced explicitly through the 
external potential. In practice, such coupling can be 
realized by scaling up or down the pseudopotentials or 
nuclear charges of initial and final molecules at their dis¬ 
tinct clamped geometries. Note that kinetic and poten¬ 
tial electron energy terms are only implicitly dependent 
on A, namely through the electron density’s dependency 
on the A-dependent external potential—which is imposed 
through application of the variational principle. 


Alchemical teleportation of an atom 

To illustrate the idea of alchemical changes within 
quantum chemistry, we now consider a process which is 
trivial when done through a realistic path, and non-trivial 
when done alchemically: The “teleportation” of an atom 
from one site to another with the constraint that the total 
number of electrons and protons is kept constant. Thus, 
instead of the trivial real space displacement of the atom, 
we continuously decrease the nuclear charge (annihila¬ 
tion) at one site while continuously increasing (creation) 
the nuclear charge at the other site by the same amount. 
For example, the external potentials of an atom at two 


sites can be linearly coupled through an alchemical path, 
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H(A) = r + r« + z£(<i^> + Ir A_), ( 6) 

where the respective atomic sites are located at the origin 
and at R. Considering only the endpoints (A = (0,1)), 
the location of the atom obviously shifted from origin 
to R. For any intermediate value of A, however, the 
electrons will distribute among the two competing poles 
of the external potential given in Eq. forming an 
attractive chemical bond. 

To numerically exemplify this process, we have chosen 
the highest occupied molecular orbital (HOMO) eigen¬ 
value, e, as property of interest, and an alchemical change 
corresponding to the linear teleportation of a Z = 2 and 
N e = 2 system, i.e. effectively translating the He atom. 
The numerical calculation of e for variable A has been 
carried out using pseudopotential interpolation within 
plane-wave basis set PBE DFT calculations, in analogy 
to previous studies. [MUCH] See Sec. | for more details. 

In Fig. [2ja), the A-dependence of £ is shown for various 
distances between the two atomic sites, d = |R| Clearly, 
while alchemical paths for small d yield simple parabolic 
shapes of s, for teleportation involving larger interatomic 
distances £ develops into a double hill. £ versus d is plot¬ 
ted in Fig. Hb) for various A values. We note that for A 
= 0.5 (magenta), the d dependency of £ corresponds to 
the case of stretching H 2 . £ increases monotonically at 
A = 0.1 and A = 0.2 as d increases. For these A values, 
the buildup of integrated electron density at the R, is still 
negligible, Fig. [2jc). Overall, the effect of nuclear poten¬ 
tial in Eq. (p|, | r 2 A R | , amounts to a static electric held, 
which induces static Stark effect. [MHlEj Because the elec¬ 
tric held decreases according to Coulomb’s law oc g, 
£ rises as a result of decreasing electric held perturba¬ 
tion. Apart from the delocalization error of DFT.[49 ) i50 ] 
such nonlinear behavior could also be related to the 
instability of H^-like systems, which has been shown 
analytically. m Hogreve pointed out that strongly po¬ 
larized electron density of asymmetric H^-like molecule 
severely destabilizes the system. [52] While the additional 
electron stabilizes the system, nonlinear behavior can be 
expected for £ in the case of strongly polarized density, 
i.e. for A > 0.3 (in Fig. He)). Fig.H c) displays integrated 
electron density slices, A(z) = f dxdy n(x,y,z), for var¬ 
ious A values at interatomic distance, d = 5A. Note that 
for A = 0.5, the electron density distribution corresponds 
to H 2 . The non-linear dependency of electron density 
n on linearly changing growth of nuclear charge can be 
seen in Fig. [2](d) for the abrupt changes in electron den¬ 
sity response induced by going from A ss 0.2 to A « 0.3. 
To investigate the impact of parameterized exchange cor¬ 
relation potentials in DFT, Cohen and Mori-Sanchez cal¬ 
culated similar changes for N e = 1 and N e = 2 using the 
hydrogen atom plus one additional atomic site where a 
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FIG. 2. Alchemical transportation of a He atom, (a) e 
as a function of A for various distances d £ {1,2,4, 7}A de¬ 
noted by solid, dashed, dash-dotted, dotted, respectively, (b) 
e as functions of d. (c) Integrated electron density, A(z) = 
J dxdy n\(x,y,z) for various A at d = 5A. The electronic 
cusps at the nuclear sites have been highlighted by their cor¬ 
responding A symbols, (d) Integrated response of electron 
density due to changing A, dA(z) = f dxdy d\n\(x,y, z) for 
various A at d = 5A. 


nuclear charge is grown, i.e. Z(A) with Z (A = 0) = 0, 
Z (A = 0.5) = 1 (H), and Z (A = 1) = 2 (He). [53] 


RATIONAL COMPOUND DESIGN 
Motivation 

The goal of rational compound design (RCD) cor¬ 
responds to solving the inverse question, i.e. “which 
compounds exhibit a set of pre-dehned desired proper¬ 
ties?” , at a rate that is superior to mere screening. [5] [SH¬ 
UT] Various approaches tackle this problem, including 
the inverse spectrum approach.[58] linear combination 
of atomic potentials,^, 6U| and many others.1561 IUlt - 
163] For the electronic potential energy, an alchemical 
path coupling E (A = 0) of one molecule to an unknown 
E (A = 1) of another compound makes explicit the com¬ 
positional dependence of the energy. Understanding such 
a dependence holds promise to dramatically reduce the 
computational burden of having to stubbornly screen one 
compound after the other. More specifically, we can ex¬ 
pand E in A in terms of a Taylor series, 

E( A) = Ei + Ad \Ei + —d\Ei + ■ • • , (7) 

where the subscript of Ei represents the quantum me¬ 
chanical expectation value of Hi. In other words, if all 
derivatives of Ei were available one could simply follow 
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a steepest descent procedure to screen a set of coupled 
“neighboring” molecules, e.g. with small differences in ge¬ 
ometry or stoichiometry, to identify and proceed to more 
promising compound candidates. Fig. 0b) illustrates the 
exploration of CCS following such alchemical predictions. 
Ideally, only a single calculation of the electronic ground- 
state Ei would be required (denoted by black circle). The 
energy of neighboring compounds can then be estimated 
via Eq. ([8]) (denoted by white circles). As we discuss be¬ 
low, it is possible to make such scans through changes in 
geometry as well as composition. 

In Ref. Q2j we already discussed that for any iso- 
electronic alchemical change, the first order derivative 
is simply the Hellmann-Feynman derivative. ]64j Conse¬ 
quently, differentiation of Eq. (§ yields, 

d\E[nx, X] = (d x H) x = J dr n x (r) <9 a ^(i\ A), (8) 

which is the same as the first order perturbation term. pH)] 
Higher order derivatives can be evaluated or approxi¬ 
mated by linear response theory. [HKII67] and will be dis¬ 
cussed below in the context of linearizing the energy in 
A in Sec. (|. 


Alchemical changes in geometry 


We now consider alchemical changes that only in¬ 
volve teleportation. To demonstrate the versatility and 
transferability of the discussed approach, we have cal¬ 
culated alchemical predictions of changes in binding en¬ 
ergy for two very different modes of binding: The cova¬ 
lent interatomic potential in hydrogen fluoride, as well 
as the hydrogen-bond-dominated van der Waals poten¬ 
tial of the water dimer. In both cases the binding energy 
is given as the difference in potential energy of dimer 
(dim) and (relaxed) monomers ml and m2, E b (d) = 
Edim{d) — E m i — E m 2 . Any approximate solution of the 
electronic Schrodinger equation at some initial distance 
di enables us to estimate the binding energy of any other 
d using the Hellmann-Feynman derivative and first or¬ 
der Taylor expansion in the alchemical teleportation path 
(Eqs. ®), 

E b {d) » E^(d) = EtW+dxEbidi). (9) 


Considering now the case of di corresponding to the 
equilibrium distance, d eq , the insets of the two top pan¬ 
els in Fig. [ 3 ] show the resulting scatter plots of Ej [ (d) 
versus the actual E b (d) for various values of d in the 
case of HF and (H 2 0)2- While there is clear correlation, 
the scale differs dramatically for the two modes of bind¬ 
ing. Most importantly, in the case of the dissociative 
tail E T1 correlates practically linearly with the actual 
binding energy. Consequently, if we now approximate 
the true E b ps E% = aij r E^ x + bi / r , (l and r corre¬ 
spond to the left-hand repulsive wall and the right-hand 



FIG. 3. Actual (black lines) and alchemical (blue squares) 
binding energy E b of repulsive (filled) and attractive (empty) 
regions of binding potentials for HF (a), (HaO )2 (b), HC1 
(c), and H 2 O-HF (d). Each screen corresponds to using only 
one self-consistent field (SCF) calculation at di = d eq , to¬ 
gether with the first order Taylor-expansion based model, 
E b = ai/ r E f 1 + bi/ r (Eq. (0). Insets in (a) and (b) show 
Eb versus E b x . The screens in (c) and (d) are slightly less 
predictive because they are made using SCF results from HF 
and (H 2 0) 2 , respectively. d eq is set to 1, 2.8, 1.4 and 2.8 A 
for (a), (b), (c), and (d) respectively. 


attractive tail, respectively) one can solve for the coef¬ 
ficients if further constraints are known. Since this is 
a rather exploratory study, we here simply assume that 
(i) E b (d = d eq ) = E^ 1 (d eq ), and (ii) E b [d -A 00 ) = 0 
in the dissociative region of the curve, and (iii) in the 
case of the repulsive region that E b (d = |d eg ) = 0 for 
covalent binding, and E b {d = | d eq ) = 0 for intermolec- 
ular binding. Assumption (iii) is based on experience 
using typical Morse and Lennard-Jones parameters. All 
resulting coefficients {ai/ r , fo;/ r } are specified in Ref. 68] . 
The predictions for scanning the entire binding potential 
agree reasonably well with the true binding potentials, 
and are shown together for both systems in the top pan¬ 
els in Fig. [3] Integrated deviations of these predictions 
are also shown in Table [I] yielding single digit percentage 
error for predicting the integral over the covalent bond¬ 
ing potential of hydrogen fluoride, and ~14% error for 
the integral over the van der Waals potential of the wa¬ 
ter dimer. We stress that the entire screen using this 
model only requires a single self-consistent field cycle to 
calculate energy and derivatives at d = d eq . 


Alchemical changes in stoichiometry 

We now extend the use of Eq. (0 to also make pre¬ 
dictions not only for teleportation changes in geometry 
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TABLE I. Numerical integrals of reference energies Eb (REF) 
and of absolute deviation of alchemical predictions E£ from 
reference energies Eb (PRE-REF) over the binding region, 
i.e. for all d where Eb < 0, and percentage thereof (%) for 
the repulsive wall predictions as well as for the attractive tail. 
Columns correspond to (a) HF, (b) (HtO^, (c) HC1, and (d) 
H 2 0-HF on display in Fig. (J3J), 


Integral [eVxA] 

(a) 

(b) 

(c) 

(d) 

REF [eV x A] (wall) 

-2.542 

-0.059 

-3.054 

-0.095 

PRE-REF [eV x A] (wall) 

0.150 

0.009 

0.656 

0.036 

% (wall) 

5.9 

15.6 

21.5 

37.6 

REF [eV x A] (tail) 

-8.199 

-0.275 

-6.017 

-0.594 

PRE-REF [eV x A] (tail) 

0.692 

0.031 

0.239 

0.072 

% (tail) 

8.4 

11.4 

4.0 

12.1 


but also for transmutational changes in stoichiometry. 
In particular, we have calculated predictions for chang¬ 
ing hydrogen fluoride into hydrogen chloride at various 
distances, as well as changing the water dimer into the 
water-hydrogen fluoride complex. Since we use pseu¬ 
dopotentials for both of these changes the total number 
of valence-electrons in our calculations does not change. 
To calculate Ej 1 according to Eq. § we have chosen d eq 
to correspond to the equilibrium distance of the target 
system, i.e. HC1 and H 2 O-HF. Again, the same assump¬ 
tions (i)-(iii) as above are used to calculate au r and bu r 
to obtain a linear approximation of the actual Eb{d) in 
E ^ 1 . Also for these changes, the resulting coefficients are 
specified in Ref. j)8]. The predicted binding curves show 
reasonable agreement with the actual numbers, as shown 
for both systems in the bottom panels in Fig. [3j Again, 
integrated and relative errors are given in Table [T| and 
show a reasonable albeit slightly worse performance than 
in the case of predicting the water dimer or the hydrogen 
fluoride. We reiterate, however, that the entire screen re¬ 
sults from only one self-consistent field cycle carried out 
to calculate energy and derivative of another molecular 
system—at the d eq of the target system. While it is also 
possible to use other d to calculate energies and deriva¬ 
tives this typically leads to less accurate predictions. We 
do not think that this constitutes a problem since knowl¬ 
edge about equilibrium distances of target structures can 
easily be obtained from literature or through inexpensive 
force-field or semi-empirical quantum chemistry calcula¬ 
tions which incur negligible computational overhead. 


Linearizing chemical space 

As we have seen above for the teleportation of the 
He atom, as well as in other studies, HU EH there are 
cases when the first order Taylor expansion of Eq. ([ 9 ]) 
does not provide satisfactory predictive power. This is 
not surprising since changes in composition correspond 
to large perturbations that typically lead to non-linear 
responses. We believe that the good performance ob¬ 


tained above for the binding curves is due to cancella¬ 
tion of higher order effects and due to the calibration of 
the linear model to the appropriate physical dissociation 
or repulsion limits. One way to systematically improve 
the predictive accuracy consists of including increasingly 
higher-order terms. Sebastiani and coworkers [571 [Ml [70] 
as well as Geerlings, De Proft and others [711175] pro¬ 
posed promising approaches in this direction. For ex¬ 
ample, akin to our discussion above, Benoit, Sebastiani 
and Parrinello investigated the performance of second or¬ 
der linear response theory for screening the potential en¬ 
ergy surface of the water dimer, and achieved very high 
predictive power. [74] How to efficiently calculate sus¬ 
ceptibility accurately and in general, however, is still a 
matter of current research. Furthermore, typically one 
observes a (sometimes dramatic) increase in computa¬ 
tional cost due to wave function-dependent susceptibili¬ 
ties, thereby defying the original motivation of RCD to 
navigate CCS without having to solve Schrodinger’s equa¬ 
tion from scratch for each and every new geometry or 
molecule. As pointed out in Ref. [T~5j . a promising alter¬ 
native route towards improving the predictive power of 
the first order derivative consists of deviating from the 
assumption that the alchemical coupling must be linear 
in A. In fact, as already mentioned above in the context 
of interpolating force-fields, [21] we are free to use any 
kind of coupling as long as we meet our endpoints, i.e. 
comply with Eq. ([!]). More specifically, if we knew the 
form of some coupling external potential v ex t (r, A) that 
induces such changes in the electron density that E{ A) 
becomes linear in A, then Eq. (|9]) would result in per¬ 
fect predictions. The quest for such a potential has been 
discussed in Ref. [12], in particular in connection to a 
1-ounce-of-gold prize for anyone who provides a solution 
to this problem. 

For a coupling path to generally fulfill the requirement 
that E( A) becomes linear in A we note that the potential 
must have such a shape that the first order derivative, 
d\E is a constant (as already pointed out and used in 
Ref. [13]), and that furthermore, all higher order energy 
derivatives must be zero. Consequently, 

0 = d™E = J dr d™~*{n\(r)dxv ext (r,\)), (10) 


V m > 1. This imposes certain constraints on the in¬ 
terpolating potential. For example, in the case of the 
second order derivative, equating the integrand to zero 
and solving for the electron density’s response results in 


d\n(r) 


~n\(r) 


d\Vext{r, A) 
d\v ext (r, A)' 


( 11 ) 


Similar expressions can be obtained for higher order den¬ 
sity response functions. Possibly, Eq. (10) could be trans¬ 
formed into a variational problem that yields an interpo¬ 
lating potential with the desired effect that the associated 
energies that are indeed linear in A. 
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CONCLUSIONS 

We discussed recent theoretical developments and ap¬ 
proaches based on coupling states using unrealistic “al¬ 
chemical” paths. Numerical evidence has been pre¬ 
sented for the applicability and versatility of alchemi¬ 
cal approaches applied to the inexpensive prediction of 
quantum mechanical observables of novel systems. The 
derivative based predictions certainly reflect the quali¬ 
tative trend of the desired binding potentials, and are 
accurate within single, or low double, digit percentage 
accuracy. Results, discussions, and current state of the 
field indicate that the study of generalized coupling ap¬ 
proaches still holds great promise for the predictive sim¬ 
ulation of molecular and materials properties, as well as 
for rational compound design. 
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COMPUTATIONAL DETAILS 

All calculations have been carried out using Kohn- 
Sham DFT[75l as implemented in CPMD[76j with PBE 
(He teletransportation) or PBEO (all other calculations) 
functional PDI [771 . Goedecker pseudopotentials [751450] 
have been used as published by Krack, [5Ij , in conjunction 
with 100 Ry plane-wave cutoffs in isolated 30 x 15 x 15 A 3 
box for He, 20 x 20 x 20 A 3 for HF—>HF and HF^HCl, 
and 110 Ry plane-wave cutoff with isolated 25 x 15 x 15 A 3 
box for H 2 O—>H 2 0 and H 2 O—»HF. Alchemical coupling 
has been imposed through linear interpolation of corre¬ 
sponding pseudopotential parameters 0 EH EH] er(A) = 
(Ji + \(<Tf — cq), where cq and 07 represent the parameters 
for atoms in Hi and Hf respectively. HOMO eigenvalues 
e have been calculated as a finite difference relying on 
Jana‘k and Koopman’s theorem,[52] [53] e a; En + s ~ En ; 
where 5 is 1% of a positive unit charge. The geometry 
scans of HF and HC1 have been performed by fixing the 
heavy atoms at origin, moving H in d direction, while 
aligning the HF or HC1 bond along d-axis. In the case of 
the (H 2 0 ) 2 scan, all geometries have been relaxed, set¬ 
ting the oxygen of the H-acceptor at the origin, while 
aligning the O-H bond of the H-donor with the d-axis. 
The H 2 O-HF geometry scans have been performed by re¬ 
placing the oxygen of the H-donor by F and annihilating 
the other hydrogen while keeping the HF bond aligned 
with the d-axis. 
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